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I. 



INTRODUCTION 



Bi-orthogonality relationships for porous media have been developed by Pro- 
fessors Scandrett and Frenzen [Ref. 1]. They showed that Fraser’s [Ref. 2] bi- 
orthogonality relationship for the propagation of Rayleigh-Lamb modes in a plate 
with traction-free surfaces is a special case of the porous case in which the porous slab 
has zero porosity. Bi-orthogonality relationships can be used to solve scattering prob- 
lems from interfaces of discontinuity between vertically stratified porous/elastic/fluid 
media. In this thesis, we will examine the bi-orthogonality relationship for the elastic- 
fluid case. 

Before the bi-orthogonality relationship can be used, the eigenvalues and eigen- 
functions corresponding to the boundary value problem must be found. In this thesis, 
we will attempt to find these numerically. The question we ask ourselves is which 
numerical approach is the best. One method is to determine Rayleigh-Lamb modes 
[Ref. 3] in the elastic solid and try to perturb those solutions to obtain modes for 
the fluid loaded layer. The problem with this method is that modes may be missed. 
Another approach is to use the Rayleigh- Ritz method [Ref. 4], which is based upon 
the calculus of variations. It uses an expansion of elementary functions to “minimize” 
a functional to the corresponding boundary value problem. The functions must be 
carefully chosen due to the complexity of the resulting coefficient matrix which will 
have terms involving the unknown eigenvalues. 

The method chosen in this thesis is the finite difference method. This method 
is very useful when numerically modeling boundary value problems involving deriva- 
tives. Using the finite difference approximations, we produce a discrete linear operator 
which approximates the differential operator of the boundary value problem. We can 
then use MATLAB [Ref. 5] to solve for eigenvalues of the linear operator which 
are theoretically close to the eigenvalues of the differential operator. We will also 
use additional numerical methods, to be discussed later, in an attempt to refine our 
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approximate eigenvalues. 

Chapter II will provide a background in acoustical theory and will introduce 
terms that the reader will need to be familiar with for further chapters. Also, the 
concept of bi-orthogonality is explained in detail and a simple example is presented. 
Finally, a description of numerical techniques and methods that are used throughout 
the thesis will be given. 

In Chapter III, background on the elastic-fluid problem is introduced. Detailed 
analytical solutions to several elastic half space problems are provided to establish the 
basis of the numerical solution. The fluid loaded plate case will be the case studied 
in the numerical portion of this thesis. 

Chapter IV examines the numerical diiTerence approximations used in the 
problem. The modal wave speeds are presented and compared to the ‘defined” wave 
speeds. Comparisons will also be made between the analytical and numerical eigen- 
functions/vectors of the fluid loaded elastic plate. 

Finally, Chapter V will provide conclusions and discuss areas of further re- 
search. 
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II. THEORY AND METHODS 



A. ACOUSTIC THEORY 



1. Basic Acoustics and the Wave Equation 

Acoustics is generally associated with sound. However, the definition of acous- 
tics is the generation, transmission, and reception of energy that is in the form of 
vibrational waves. Audible waves are between 20 and 20,000 Hz where \Hz = 
\ cycle /second. The simple oscillator [Ref. 6], or basic spring mass problem, is the 
most common example used to demonstrate basic acoustics. The simple oscillator is 
governed by the equation 




(III) 



with a general solution of 



X = Acos{uiot) + Bsin{u}oi) 



(II.2) 



where cjo is the angular frequency in radians per second (rads/sec). Since there are 
27 t radians per cycle, we define frequency as 




(II.3) 



In addition, we note that the period T is equal to -j. If the initial velocity and 
displacement are known, we can obtain an exact solution to Equation (II.l). 



The fluid-elastic problem that is the central problem of this thesis is a two- 
dimensional problem, and thus the fluid and elastic potential functions are governed 
by the two-dimensional wave equation 




(II.4) 
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We assume an incident wave is travelling in the positive ^-direction, which strikes an 
interface generating reflected and transmitted waves. A plane time harmonic wave 
propagating in the postive x-direction has the form 



$ = (II.5) 

where k\ = — . Substituting Equation(II.5) into Equation(II.4), one finds the ampli- 
tude function must satisfy the relation 



+ (^'i - = 0 (n.6) 

which is the canonical form of the boundary value problem to be solved. 

2 . Stress-Strain Relationships For Elastic Solids 

Elastic solids have the property that when external pressures are removed, the 
solid returns to its natural state. Strain is a measure of displacement while stress is 
a measure of force per area. We define the strain as 



1/ N 

Sij — 



(II.7) 



where u is a displacement vector and tensor notation has been adopted. The term 

Uij means the partial derivative with respect to Xj, or for example, u \^2 corresponds 

to — . The strain, Sij^ corresponds to elements of a matrix. Note that Sij is dimen- 
0x2 

sionless. Stress, has the dimensions of force per area. For small strain, 5 , the 
stress, r, is proportional to strain. This relationship is known as Hooke’s Law [Ref. 
7], which can be written for an isotropic homogeneous solid as 



Tjj — ^^kk^ij T (1 1. 8) 

where A and /i, known as Lame’s elastic constants, have the dimension of pressure, 
and Sij is the well-known Kronecker delta defined as follows: 
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( 11 . 9 ) 



Sij — 



1 t = J 

0 i^j 

The displacement vectors can be represented in terms of potential functions. 
The displacement equation of motion for an isotropic homogeneous elastic solid [Ref. 
7] can be written as 



+ (A + /li) VV • u = pii 

We consider a decomposition of the displacement vector in the form 



(II.IO) 



u = V(^ + VxV’ 

It can be shown that the displacement vector u satisfies Equation (II.IO) if 



1 - 

VV = - 2 ^ 
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( 11 . 12 ) 



where ^ and = —[Ref. 7]. The question of the completeness of 

, P , P 

the solution in Equation (II. 11) is guaranteed by the following theorem given in 
Achenbach’s book on elastic solids[Ref. 7]. 

Completeness Theorem: Let 



u € C^V X T) 
f € C(V X T) 



and satisfy the equation 



+ (A + /x) VV • u + pf = pii 

in a region of space V and in a closed time interval T. Also let 
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f = 4 vF + 4V X G. 



Then there exists a scalar function and a vector-valued function in 

terms of which u(x, t) is represented by 



with 



u = V X 



V • V’ = 0, 

and where v^(x, t) and rf){x,t) satisfy the inhomogeneous wave equations 

VV + 

+ G = 

Op 

The proof of this theorem can be found in [Ref. 7]. 

Representation of the displacement vector u in terms of potentials will become 
important in the remainder of this thesis. We restrict our analysis by assuming the 
strain in the z direction is zero, which corresponds to the two dimensional plane strain 
case. This allows us to specialize our decomposition into the following 

U = (Px + i’y 
V = (Py-i)x 

where the vector valued displacement potential is reduced to a single nonzero scalar 
component. Our stress equation follows from Equation(II.8) as 

Tj;y = jJ.{Uy + V^) 

Tyy = -f Uj,) + 2jJ,Vy 

where, as before, the potentials (p and ip satisfy Equations (11.12) 



(11.14) 



(IU3) 
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For an ideal fluid the stress in any direction equals the pressure stress of the 
fluid, or in equation form 



Trx — '^yy — P 



(11.15) 



B. BI-ORTHOGONALITY 

We say that two functions t^(x) and V>(x) are orthogonal over an interval 

0 < X < L, if / c^(x)i/i(x)dx = 0. This is analogous to the dot product of two 
JQ 

orthogonal vectors equaling zero . The orthogonality relationships[Ref. 8] for sines 
and cosines are as follows 



I 



. riTTX . mTTX 
sin sin — - — ax = < 

0 L L 



0 n ^ m 



^ n = m 



i 



mrx mTTX 



cos—— 
0 L 



cos- 



dx = ^ 



0 n ^ m 

^ n = m 7 ^ 0 

L n = m = 0 



i 



^ . nirx mirx 
sin — —cos — - — dx = 0 



(11.16) 



(11.17) 



(11.18) 



10 L L 

This orthogonality relationship allows us to determine values of coefficients when 
arbitrary functions are expanded in terms of these functions. For example, if /(x) 
can be written in the form 



/(^) = Bnsin 



nwx 



(11.19) 



71 =] 



then we can apply the orthogonality of sines over the interval [0, L] to obtain 



/ f(x)sin'^dx = Y:Bj si 
Jo L Jo 



nnx . mxx 



sin — —sin- 



dx 



( 11 . 20 ) 



which leads to 
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or the more familiar form 



Bn = ^ 



nnx 
sin — —ax 



I 



sin 



<ynr^x , 

^ dx 



L 



( 11 . 21 ) 



2 nTTx 

Bn = — f(x)sin — —dx (11.22) 

L Jo L 

The analysis is similar for cosines. Note that we have formally interchanged summa- 
tion and integration in the above steps which is valid provided the series and integrals 
converge. 

Bi-orthogonality relationships generalize the orthogonal relationships of clas- 
sical eigenfunction expansions. In [Ref. 1], these relationships were derived for waves 
propagating at a free surface and along interfaces between porous and elastic media. 
The general bi-orthogonality relationship for a two dimensional porous layer which is 
fluid loaded is: 



/ ~ + s^U^].porousdy + / \-p^U^]siuiddy = 0 (11.23) 

where the fluid layer and porous layer thicknesses are h\ and h 2 respectively, and 
s — —[dp where f3 denotes the porosity and p denotes the interstitial fluid pressure. 
Note also that A and B denote different modes within the same solid. When the 
porosity of the elastic media is zero, the equation reduces to the fluid loaded elastic 
solid bhorthogonality relationship. 

To illustrate how the bi-orthogonality property can be applied, we will examine 
the simple case of two fluids, each with different densities in contact at an interface 
along the y axis. In this case, the bi-orthogonality relationship reduces to standard 
orthogonality between sines and cosines. A diagram is shown in Figure(l), 
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y 


II 

O 


II 

O 


Fluid 1 


Fluid 2 


In 


Tn 


Rn 




Py = 0 


0 p = 0 



Figure 1. Bi-orthogonality of Two Fluid Layers 

where fluid 1 satisfies homogeneous Neumann[Ref. 8] conditions, py = 0, and fluid 2 
satisfies homogeneous Dirichlet conditions, p = 0, along the lateral sides of the layers. 
From [Ref. 8] we know the eigenfunctions are sines for the Dirichlet conditions and 
cosines for the Neumann conditions. We equate the pressure and flow normal to the 
interface to obtain 



n 1 




{In + Rn)p^n\y) 




where p^^{y) and U!^'^{y) are defined as follows; 



Pn\y) = cos(ny) and p^^{y) = sin(my) 
U^\y) = cos(nt/) and U!^\y) = sin(my) 
Applying the bi-orthogonality condition as follows 



(11.25) 



uP\p' = p‘]dy 

£ pP\U' = U^]dy 

j( t'fV' = p‘]dy 

r = U^]dy, 

Jo 



(11.26) 
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we obtain the following two sets of equations, written in matrix form for easier ma- 
nipulation 



J(/-h^) = KT, 

J(7 - R) = G^f, 

( 11 . 27 ) 

G{I + R) = QT, 

K^(7 - R) = Qf. 

We now define Q = (Q~*G) and K. = (J“’K), which, if both formulations are to give 
identical answers, results in an identity = I. Solving for ^ and /C we have: 



/ co${py)sin[ny)dy 

Qrn = (Q-'G) = - —f — 

J sin^(ny)dy 

f cos{qy)sin{ny)dy 
fCn, = (J-'K) = ^^2—^ 

/ cos^{qy)dy 
Jo 



leading to 



( 11 . 28 ) 



Gpn — 



n 



7T \n^ — {p — ly ^ 



where p + n must be even 



2(2-S„) f 



where q + n must be even, 



( 11 . 29 ) 



which gives the matrix element as 



^ ^ Gpn^nq — 



0 if p + 9 is odd 
S{2-S,,) 






n 



n - (P- - (?- 1)^] 

It can be shown that the above summation gives 



if p -t- 5 is even. 

( 11 . 30 ) 
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{g)CU = 6^, (11.31) 

Knowing these relationships and given initial incident wave amplitudes, /, we can 
solve for the reflected and transmitted amplitudes from either one of the pairs of 



vector equations given by Equations (11.27). 



C. NUMERICAL METHODS 
Newton’s Method 

Newton’s method, one of the numerical methods used in this thesis, is a widely 
used method for solving nonlinear equations. We will briefly review the concepts and 
give a short example of this method. Newton’s Method uses a line tangent to the 
curve of a function and is based on the linear approximation of the function. We start 
with an initial guess Xq, that we believe to be close to the root, and determine f{xo). 
Then we move along the tangent line until the point intersects the x-axis. This value 
of X now becomes our next iterate. In general terms, Newton’s Method is given as 
[Ref. 4] 



f{Xn) 

^n+1 — rf/ 

fi^n) 



(11.32) 



where n = 0, 1, 2, — Newton’s Method converges rapidly when the initial guess is in 



the neighborhood of the root. In fact, Newton’s Method is quadratically convergent. 
Here is an example [Ref. 4]; 



f{x) = 3x + sin X — e® 
/'(x) = 3 + cos X — e® 



We choose xq = 0 and perform 3 iterations as follows 

xi=xo- 4^ = -33333 

X2 = o:, -4^ = .36017 

X3 = X2- 4^ = .3604217 

f{^2) 
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which is accurate to 7 significant digits. Although the convergence is rapid, there are 
other considerations that need to be addressed with this method. Newton’s Method 
requires 2 function evaluations per step. However, the total number of function 
evaluations is generally the same or less than other methods and Newton’s Method 
usually takes fewer iterations to determine a root with high accuracy. The main 
drawback to Newton’s Method is that it will not converge for certain functions. This 
happens when the tangent line is approximately horizontal or when we repeat a value 
previously used, for example Xi = X 4 . Although these are strong considerations, 
this happens infrequently. For our needs, we will use this method to refine our wave 
speeds found in Chapter III. We will also use this method in refining our eigenvalue 
estimates found using the finite difference method which is discussed in Chapter IV. 

2 . Finite Difference Method 

The finite difference method is used in numerically approximating derivatives 
and integrals of functions. We can use this ability to replace derivatives with difference 
quotients and thereby approximate solutions to difficult differential equations. There 
are three basic types of difference quotients; backward, central, and forward, with 
central being the most accurate. In certain instances, such as endpoints, backward 
or forward differences are often used due to the lack of additional points beyond 
the extremity of an interval necessary for accomplishing the central difference. In 
any case, it is customary to ensure that all difference operators have the same or like 
order of accuracy to avoid extra work and increase the possibility of error cancellation. 
Difference operators can also be applied more than once in order to obtain second 
order and higher difference approximations to higher order derivatives. For this thesis, 
we will concern ourselves with only second order differences. The central difference 
operators [Ref. 4] are as follows: 
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(11.33) 







2h 



+ 0{h^) 



d X. 2xj 4“ 1 /o/l2\ 

= p + 0(A^) 

where /i = Aa; is the stepsize. Since forward and backward differences are 0{h), we 
need to expand Xi+i and x,_] as follows to obtain formulas of higher order for our 
backward and forward difference operators: 



Xi+i = Xi + x\h + x"/i^ + x'l'h^ H 

x,_] = Xi — x'-h + x"K^ — x'l'h^ + • • • 
Re-arranging terms we obtain 



(11.34) 



X 



/ 

i 



X 



/ 

i 



Xj'+i X{ 




h 



X • 

2 ^ 

-b -x'! 
^2 * 



- yx” + 0(A») 

-^xf+0(A») 



forward difference 
backward difference 



(11.35) 



We will be able to substitute for the terms x", and x'" in our final discretization. This 
will be discussed in more detail in Chapter IV. 
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III. DISPERSION EQUATIONS FOR 
ELASTIC HALF SPACES 



In this chapter we will examine several elastic fluid problems with infinite 
and finite layers. In each problem a relationship known as the Dispersion Equation 
is found. A system is dispersive if the phase velocity of a wave depends upon the 
frequency of the wave. Solving this dispersion relationship at a prescribed frequency, 
we can determine the wavespeeds for propagating modes and use this for an analysis 
of the system. These modes are also the necessary constituents in applying the bi- 
orthogonality relationships being studied. 

A. INFINITE ELASTIC HALF SPACE WITH A FREE 
SURFACE 



Free Space 



Infinite Elastic 
Half Space 



Figure 2. Infinite Free-Elastic Half Space 

The case of an infinite elastic half space with a free surface is described by 
Achenbach [Ref. 7]. We consider the two dimensional plane strain case. We apply 
stress free boundary conditions at the free surface: 

jX^Uy -j- 0 [Ill.l^ 

'^yy — + ^3^) + = 0 (HI. 2) 
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Dividing Equation (III. 2) by p and substituting our longitudinal and transverse wave 
speeds, and Cj we obtain 



Assuming a solution for our vector potentials, (p and -0, as follows, 
our solutions for u and v become 



(III.3) 



(III.4) 



U = + Xl)y = ikcp + ijjy 

V = P’y — tkx = — iklj} 

and the two boundary conditions can be written: 



(III.5) 



Uy + ikv = 0 



(c^ — 2cj)iku + cj^Vy = 0 
Next, our potential functions must solve Helmholtz equations in the elastic solid 



(III.6) 



(II1.7) 



V V +hly? = 0 

Dropping the tilde notation on the “amplitude”functions of y and the 
pendence, these functions satisfy the wave equation if 



= 0. 

0"-f0 = O, 

where = k^ — kj^ and — k^ — k^. This leads to exponential solutions 



(III.8) 



= Ae~’”‘ + Be™ 
0 = Ce~^^ + 



(III.9) 
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where we stipulate that Re{^, rj) > 0. To ensure boundedness of the solutions as 
y — ^ — oo, we neglect solutions that “blow up” to obtain 



ip = 

rp = 



Substituting into the boundary conditions we have 



(III.IO) 



ikrjtp + + ik{rjip — iktp) = 0 

{cl - 2(^)ik{iktp + + cl{i]'^(p - ik^xp) = 0 

Evaluating the amplitude functions at ?/ = 0 and simplifying leads to 



(III.ll) 



{2k^ - kl)B - 2ikiD = 0 
2ikriB + (2P - kl)D = 0 



( 111 . 12 ) 



The determinant of the coeflSicient matrix in Equation (III. 12) leads to the Dispersion 
Equation 



(2^2 - klf - ^k'^iy = 0 (III.13) 

In this instance, the waves are nondispersive because the phase speed determined 

by the above relationship is independent of the frequency. The frequency can be 

scaled out by noting k = —, and dividing the equation by Using cl = 5690m/s 

c 

and ct = Zl4:5m/s, the compressional and shear wave speeds for steel, we can apply 
Newton’s Method to solve for Cr, the Rayleigh wave speed, and find 

cr = 2907mls 

where cr < ct < cl. From [Ref. 7] we have the following approximate formula for 
computing the Rayleigh wave speed. 



CR 



.862 + 1.141/ 
1 + 2 / 



-Ct 



(III.14) 
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Using this formula, we predict a Rayleigh wave speed of 2902.25m/s which is very 
close to our result. Now letting Z) = 1 we find that B = .6748U At the frequency 
u> = \ and setting x = y = 0, the real parts of the particle displacements (u,v) 
have been plotted in Figure(3). Note that this is the time history of the particle 
displacement at the surface of the elastic half space. In addition, as depth increases, 
the amplitude of the Rayleigh wave decreases and the motion alternates between 
retrograde and prograde. Rayleigh waves for various elastic solids appear in Table (I) 
later in this chapter. 




Figure 3. Rayleigh Wave Propagation 



B. FLUID LOADED INFINITE ELASTIC HALF SPACE 

Infinite Fluid 
Half Space 

I „x 

Infinite Elastic 
Half Space 



Figure 4. Infinite Fluid-Elastic Half Space 

When a fluid is added, we have slightly different boundary conditions. Equa- 
tion (III.l) still holds, but Equation (III. 2) is modified due to the fluid pressure at 
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the boundary. This leads to 



Tyy = \{Ua: + Vy) + 2fiVy = -p (III.15) 

for which 

p = -ps<i>t (III.16) 

and where ^ is the acoustic velocity potential for the fluid and p/ is the fluid density. 
The normal stress boundary condition becomes 

(?jVy + (c? — 2(3r)ux + ioj—(j> = 0 (111.17) 

Ps 

We allow slippage along the fluid/solid interface (y = 0), but demand continuity of 
the vertical velocity there, i.e. Vgoiid = This leads to a third boundary condition 

— — LOkxl) + = 0 (III. 18) 

The velocity potential function must also satisfy a Helmholtz equation 

V2<^ + 4<^ = 0 (III.19) 

UJ , . 

where kp = — and cp is the fluid wave speed. As before we determine the form of 
cp 

our solution using the wave equation and assuming a potential of the form 



<f> = (III.20) 

Dropping again the tilde notation and the dependence, the amplitude func- 

tion must satisfy 



<i>yy — 4> = 0 



(III.21) 



where = k^ — kp. The amplitude function has the exponential solution 
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<f> = Ee~^y + (III.22) 

Disregarding solutions that “blow up” as y ^ ±oo we obtain the amplitude functions 
for the solid and fluid 



(p = Be^y 

^ = De^y (III.23) 

(j) = Ee~^y 

Substituting these into the three boundary conditions we have at y = 0: 



2ikr)<p + + k^)ip = 0 

{cly^ — — 2c^))ip — 2ik^(^ij} + iu}^(j> = 0 

iu'qp + ukxl^ — (^(f) = 0 

Simplifying we obtain 



(III.24) 



2ikrjB + {2k^ — kj)D = 0 

(2fc^Cj — k\(?i^)B — 2ik^c^D + ii^^E = 0 (III. 25) 

iurjB + ukD — (E = 0 

Computing the determinant of the matrix and simplifying results in the Dispersion 
Equation. However, unlike the previous case, we now have an added term due to the 
presence of the fluid half space. This equation is 



[(2P - k^f - = 0 (III.26) 

C Ps 

Note that once again the frequency can be factored out of the above expression which 
implies a nondispersive media. Applying Newton’s Method using the wave speeds for 
steel we determine what is referred to as the Scholte[Ref. 7] wave speed, 

c, = 2903.567 - 35.325f m/s 
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Note that in the case with the fluid half space we have a complex wave speed. Letting 
E = 1 we find B = 2.51 + .395i and D = —.030 + 3.89L Setting cu = 1 and plotting 
the real part of the horizontal and vertical displacements at x = y = 0, we obtain the 
particle displacement history as this Scholte wave passes through the point x = y = 0. 
The plot looks the same as Figure (3). 

C. FINITE FLUID INFINITE ELASTIC HALF SPACE 



Fluid 

L___I 

Infinite Elastic 
Half Space 



Figure 5. Finite Fluid Infinite Elastic Half Space 

We next consider a fluid layer of finite thickness overlying an elastic half space. 
The boundary conditions from before, Equations (III. 15), (III. 17), and (III.18) still 
hold. However, an additional boundary condition is needed at the boundary between 
the fluid and free space which is that the pressure be zero there. From Equation 
(III. 16), this leads to 



p = i(jjpf(j) = Q ^ (f> = 0 (III. 27) 

Also, from before we obtain the solutions to the Helmholtz equation where 

ip — 
tp = 

Our third potential changes because we now have a finite layer 
discard one of the solutions, but instead must have 



(III.28) 
. We are unable to 
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<i> = + De~^^ 

Simplifying this by combining the terms, we choose 



(III.29) 



= C'sinh(C(/i -y)) (III.30) 

Notice that by making this choice we have incorporated our additional boundary 
condition so that when y = h, (f> = 0. Substituting into the boundary conditions at 
y = 0, we have 



or 



2ikri(f + = 0 

— 2(9p))ip — 2ik^c^xj} + iu^4> = 0 (III. 31) 

iojTup + Lok-tl^ — (C cosh(C(/i — y)) = 0 



2ikT]A + {2k^ + kj)B = 0 

(2fc^c^ — cj^)A — 2ik^c^B + sinh((h)C = 0 (III. 32) 

iu>7fA + cjkB — ( cosh((h)C = 0 

Computing the determinant and simplifying leads us again to a Dispersion Equation, 
which is 



[(2^2 - 4)2 - 4k^r]^] + 4^^ tanh(C/i) = 0 (III.33) 

C Ps 

Unlike the previous case, the frequency can no longer be factored out. This dispersion 
relation signals that we now have a dispersive media which is intermediate between 
the other two limiting cases. Notice that as a; — » 0, the “fluid term” is 0{u>^) while the 
solid terms are 0(a)‘‘), and therefore the wave acts like a Rayleigh wave. As a; ^ oo, 
tanh(C/j) — > 1 we obtain the Scholte limit. This makes sense physically because as 
cj ^ 0, A ~ 0( — ) becomes large and the fluid layer will appear to become extremely 

(jJ 

“thin” relative to the propagating wavelength. On the other hand if cj — »• oo, the 
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fluid will appear to be “thicker” in terms of the number of wavelengths that can fit 
onto the fluid layer, and we get the other limit of an infinite fluid layer. As in the 
previous problems we solve for the Scholte interface wave speed. A height of Im is 
used for the fluid layer in the computation. The wave speed for steel is 

c = 2906 - .001 h' m/s 

Letting C = 1 we find that A = .00947 — .104xl0~®z and B = .147x10“® — .0057h’. 
We now set cu = 1, /i = 1, and x = y = 0. This plot is the same as our Rayleigh wave 
plot in Figure (3), which was expected. Listed below in Table (I) are Rayleigh and 
Scholte wave speeds for various elastic solids. 



Solid 


CL 


Ct 


cr 


finite fluid 


finite fluid 


Steel 


5690 


3145 


2907 


2904 - 35f 


2906 - .OOh' 


Iron (Cast) 


4175 


2308 


2133 


2127 - 45f 


2133 - .0009^ 


Aluminum 


6242 


3144 


2930 


2904 - 112z 


2930 - .003? 


Glass (Pyrex) 


5637 


3297 


3026 


2993 - 153i 


3026 - .004? 


Rubber (Hard) 


2117 


864 


814 


694 


814 



Table I. Rayleigh and Scholte Wave Speeds for Various Elastic Solids 



D. FLUID LOADED ELASTIC PLATE 




Figure 6. Fluid Loaded Elastic Plate 

We consider an elastic layer of finite thickness. We examine the relationhip 
of a finite fluid media in contact with a finite elastic medium or plate. As shown 
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in Figure (6), the origin is assumed to go through the middle of the elastic solid. 
This will help later to show the symmetric and anti-symmetric modes as discussed by 
Achenbach [Ref. 7] in the unloaded, ‘‘free plate” problem. Our boundary conditions 
are slightly different with this case. We choose a thickness h\ for the elastic layer and 
a thickness /12 for the fluid layer. At j/ = — we have 



At y have 



Try — ^ 
'^yy “ 0 



(III.34) 



h. 



Try — 



^yy = -P 

solid ~ fluid 



(IIL35) 



Finally, at y = /^2 + we want the pressure to again be equal to zero. We again 
assume solutions to the wave equation. However, we now have a finite fluid and elastic 
layer, and therefore must retain all of the exponential terms. We write the solutions 
in terms of cosh and sinh in order to better see the symmetric and anti-symmetric 
portions of the solutions. We have, 



ip = A cosh(r/y) + B sinh(r/y) 
xj; = C cosh(<^y) -f D sinh(^y) 
(f) = Esinh{C{h - y)) 



(III.36) 



hy 



where h = — + h 2 , and we have taken into account the boundary condition that 
p = 0 a.t h. 

Substituting into our boundary conditions and letting p = —, v/e have 

Ps 
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(III.37) 



+ (2P _ = 0 

2ikr^<^y{^) + (2P - 4)V’(^) = 0 
{2k^-kl)^{-^)-2ikii,y{-^) = ^ 

(2P - kl)^{^) - 2ikii,y{^) + '^<i>{h,) = 0 
iu‘qipy{^) + ujkxJ;{^) - <j)y{h2) = 0 
Trying to find and simplify the determinant is difficult. We try a simpler approach by 
adding and subtracting like equations to obtain our Dispersion Equation. Adding and 
subtracting the first two equations and doing the same with the next two equations, 
we have the following 



2ikr]B cosh{ri^) + (2P — k^)C cosh(^-^) = 0 
2iA; 7/A sinh(77-^) + (2P — kj)D smh{^^) = 0 
(2P _ A;^)Acosh(7y^) -2fA:^Dcosh(^^) + ^Esinh(C/i2) = 0 (HI.38) 

(2P - 4)Bsinh(7/^) -2i/t^Csinh(^^) + ^Esinh(C/i2) = 0 

i<^VVyi^) + - (f>y{h2) = 0 

From the first two equations, we solve for A and B in terms of C and D. Substituting 
into the next two equations we obtain a relationship of C and D and find E in terms 
of C. This is not as difficult as it seems, but does require accurate bookkeeping to 
keep up with the algebra. Substituting all of these into the last equation, and after 
much algebra, we obtain the Dispersion Equation for a fluid loaded plate. 



[(2P — k^Y — 4k^T]( tanh(?7^) coth(^^)][(2P — k^Y ~ coth(?;^) tanh(if^)] 

= -P/^r^tanh(Cfe2) pp _ 

PsC 



(III.39) 

Because the dispersion relationships will involve infinitely many complex roots, due 
to the finite nature of the elastic/fluid system, Newton’s method may indeed be able 
to find one or several of the roots, but will not be efficient in finding the first N roots 
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(ordered by magnitude), simply because required starting values are not available. 
One lacks the confidence that all of the eigenvalues will be found with Newton’s 
method without sufficiently close starting values. 
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IV. 



RESEARCH RESULTS 



A. DISCRETIZATION PROCESS 

We begin the numerical analysis of the fluid loaded plate by applying finite 
difference approximations to the boundary value problem. We know that in the fluid 
and solid, the respective reduced Helmholtz equations must be satisfied. We define 
{N l)hf as the thickness of the fluid layer and (M + l)/ij as the thickness for the 
elastic layer, where (A'^ + 1) is the number of discrete points in the fluid layer and 
{M + 1) as the number of grid points in the elastic layer, and hj and are the 
stepsize for the fluid and elastic layers, respectively. In order to obtain our unknown 
modal wave number, k^, on the right hand side of the system of equations, several 
substitutions are necessary. We define tp* = ipn and also substitute 4>t = 4>- Therefore 
in the final analysis we will be solving for normal modes of the functions and 

rp*. Dropping the superscripts and subscripts to avoid clutter in the discretization, 
we have 



cP'! k)<P, = kl<Pi i = \..N-\ 

+ kl'Pi = KVi * = - 1 (IV. 1) 

ipi + k^xpi = k^xpi i = 1..M — 1 

The additional boundary conditions on the fluid loaded elastic plate from Chapter 
III, can be reduced after completing the above substitutions and applying the finite 
difference method, to the following: 



<pN 

2ip" + {2k1 — kj)if — 2xp' = ^(p 
—2k'p(p' kjtp + 2ip" = 0 

2(p" + {2kl - k^)<f -2xP' = 0 
—2kp^if' + k^ip + 20 " = 0 
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where i = 0 is the interface between the fluid and the solid. In order to incorporate 
these additional boundary conditions into Equations (IV. 1), we must solve for what 
we will call, pseudo nodes, A pseudo node is a point outside of the fluid or elastic 
solid that will be used in the differencing of a potential. This can be seen below in 
Figure (7), 




Figure 7. Discretization of Fluid Loaded Elastic Plate 

where is the pseudo node for the fluid and 9_i, S'lid are the 

pseudo nodes for the elastic solid. 

The equations for the pseudo nodes require some careful algebra and can be 
reduced to the following: 



<j>i — (j)-\ + 4’o + xl>o — [1 + (^T ~ ^i)]y*o — h,hfk\w^ipi — y>2 

‘P-i = [1 + ” ^l)]‘Po + (1 - -Ip2 + (2/ii + 

V--1 = (4 + k^h^,)rl>o - aV'i + (ft.4 - ^)<Po + (^ - 2klhs)^i - + ^4>o 



<PM+1 = [1 - - k^)]<PM + (1 - h^kl)(pM-l - <PM-2 + 2/liV'M-l + - ^hs)tpM 

'pM+i = (4 - h^k^)tpM - 3V'm-i + (t^ “ + ^^<PM-2 



(IV.3) 



These pseudo nodes are substituted into our Helmholtz equations from (IV. 1) to 
obtain equations for the boundaries t = 0 and i = M. Our discretized boundary 
equations are now entered into a MATLAB code given in the Appendix. Note that 
for the MATLAB code we must start our index at 1 as MATLAB does not have 0 
indices. A shift is required so that the first N places correspond to (j>, the next M + 1 
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places correspond to (^, and lastly, the next A/ + 1 places correspond to 0. This allows 
us to keep track of each of the three potentials. MATLAB computes our eigenvalues 
and orders them in ascending magnitude. The code also calculates the normal and 
shear stresses for use in applying the bi-orthogonality condition. Finally, various 
graphs involving our discretized solutions will be plotted against analytical solutions. 
The analytical solutions for stresses, pressures, and displacements are determined 
using the coefficients listed in the last section of Chapter III. 

B. NUMERICAL RESULTS 
1. Modal Wave Numbers 

The elastic solid used for the following analysis will be steel. We begin with 
an examination of the eigenvalues of our discretized matrix. In the discretization, 
u = 27T, N = 10, and M = 50, leading to a square matrix of size 112 from which 112 
eigenvalues are determined. Since we are solving for we must take the square root 
of the eigenvalues and place them in order. Determination of the accuracy of these 
wave numbers will involve substitution into Equation (III. 39), our fluid loaded elastic 
plate dispersion relation. A plot of the norm of the eigenvalues is shown in Figure 
(8) on the following page. Aside from the last few eigenvalues, we get a nice ordering 
of the eigenvalues. 

We attempt to refine our approximate wave numbers using several iterations of 
MAPLE coded Newton’s Method with the eigenvalues as initial iterates. The results 
are shown in Table (II). For some of our “refined” wave speeds, Newton’s Method 
attempts to converge to a root. Newton’s Method diverged for the roots marked with 
an asterik. After refinement, several of the roots show improvement in satisfying this 
equation, i.e. the magnitude of the residual of the relationship is much closer to zero. 
The fundamental root has a residual magnitude change from 10“^^ to 10~^°. The 
second and third roots change almost as much. However, the rest of the roots exhibit 
less of a change in magnitude, indicating that the starting values for the roots may 
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modes 



Norm of Eigenvalues 
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not be close enough to the actual roots for Newton’s Method to converge. 



n 


k 


Refined kn 


relative error 


1 


.001291 - .0003723f 


.001218 


0.2824 


2 


.001291 + .0003723f 


.001218 


0.2824 


3 


.008880 


.004932 


.4445 


4 


.006711 - .008003?* 






5 


.006711 + .008003?'* 






6 


.004703 - .01666?'* 






7 


.004703 + .01 666?'* 






8 


.002550 - .0318? 


.01127 - .0210?' 


.4355 


9 


.002550 + .0318? 


.01127 + .0210? 


.4355 


10 


.001664 - .0476?' 


.0007991 - .0398?' 


.1646 



Table II. Calculated versus Refined Modal Wave Speeds 



Table (II) shows the first three roots converging to a real root. The other “refined” 
kn's do not appear to converge to any of the kn's used as initial starting values. It 
may be necessary to search more than the first ten modes to help determine if a shift 
of the “refined” values to a corresponding is possible. For this reason, we will 
concentrate on the first few modes. 

2. Stress Potentials 

The following graphs show the analytical versus the discretized solutions for 
the fluid and solid potentials ip, and tp. Figure (9), displays the best results. The 
analytical solution of <f> is almost identical with the discretized solution. Note that 
the discretization code has (f>i as the first interior point and not the boundary. This 
accounts for the slight mismatch at the fluid-free space boundary. Figure (10) shows 
a simliar pattern where our solution is nearly identical, but diverges slightly as we 
near the fluid-elastic boundary. Figure (11) diverges at both ends but is on a scale of 
10“®. When viewed on a small scale, the discretized solution does diverge more than 
4> and (p. It is obvious that there is a problem at the fluid-elastic boundary and from 
Figure (11), there is a slight indication of a problem at the bottom boundary. 
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Analytical vs Discretized for Phi (fluid) 




Figure 9. Analytical versus Discretized Solutions for (j> 
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Analytical vs Discretized for Varphi (solid) 




Figure 10. Analytical versus Discretized Solutions for (f 
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Analytical vs Discretized for Psi (solid) 




X 10 ^ 



Figure 11. Analytical versus Discretized Solutions for -0 
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Figure (12) illustrates a pressure-stress graph. Notice that at the boundary 
of the fluid and elastic layer a discontinuity exists. One of our boundary conditions 
from the fluid loaded plate is 



T„ = -r (IV.4) 

This means that our plot should connect smoothly at the interface. The discontinuity 
indicates a problem with the discretization at the interface. This problem was evident 
in the previous plots and is re-enforced with this graph. 

Figure (13) plots the wave numbers of the first three modes as a function of 
increasing frequency. We see that our fundamental mode is purely real and our second 
and third modes are complex. 

3. Bi-orthogonality 

The bi-orthogonality condition of differing modes from Chapter II indicates 
that we should expect a diagonal matrix in forming the bi-orthogonality product of a 
set of modes, where the diagonal element, is the “norm” of each individual mode. 

= r - t' \p^'^^U^%iuudy (IV.5) 

J-h2 Jo 

The bi-orthogonality condition states that the off-diagonal elements, the product of 
two separate modes, should go to zero. Numerically, we expect these elements to tend 
toward zero. The numerical bi-orthogonality relationship for the fluid loaded plate 
yields the following matrix; 



4.441 


X 


10® 


4.561 


X 


10® 


5.233 


X 


10® 


3.639 


X 


10® 


3.836 


X 


10® 


4.561 


X 


10® 


4.441 


X 


10® 


5.233 


X 


10® 


3.836 


X 


10® 


3.639 


X 


10® 


5.118 


X 


10® 


5.118 


X 


10® 


7.046 


X 


10® 


2.670 


X 


10® 


2.670 


X 


10® 


2.617 


X 


10® 


1.005 


X 


10^ 


6.139 


X 


10® 


1.524 


X 


10" 


1.832 


X 


10" 


1.005 


X 


10^ 


2.617 


X 


10® 


6.139 


X 


10® 


1.832 


X 


10" 


1.524 


X 


10" 
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height of fluid-elastic layers 



Pressure - Normal Stress Mode 1 
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Figure 12. Pressure - Normal Stress in Layer 
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Frequency (Hz) 



Wave Propagation for Steel 
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Notice in this matrix, most of the off*diagonal elements tend to exhibit some form 
of reduction. However, this did not happen as quickly as believed. This indicates 
that the bi-orthogonality condition is very sensitive to small errors in calculating the 
modes. A more careful numerical treatment might improve the diagonal form of the 
matrix, but it may be necessary to limit the application of the numerical method 
to eigenmode determination only, from which the eigenmodes can be refined and 
ultimately found by analytical means. 
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V. CONCLUSIONS AND FUTURE 

RESEARCH 

We can draw several conclusions on the discretization of the fluid loaded elastic 
plate. The eigenvalues, or modal wave numbers, that resulted from the discretization 
nearly satisfy the Dispersion Equation (III. 39). Using Newton’s Method, we deter- 
mined that our initial values may not be close enough in some instances to the actual 
roots as is evidenced from the lack of improvement of the residual calculation in the 
“refined” root. Therefore, we conclude that there are numerous other roots that we 
have not found. 

The biggest disparity lies in the plot of pressure and normal stress versus the 
height of the fluid and elastic layer. As seen in the discretization code in the Ap- 
pendix, differencing of the potentials is necessary when solving for the stresses and 
displacements. Therefore, any numerical error in the fluid-elastic boundary differenc- 
ing is compounded. Recall, in theory the normal stress is equal to minus the pressure 
at this boundary. The discretization does not show this. Finally, we conclude that 
our graph of the first three modes in Figure (13), could indeed be modes for this 
problem. However, they are most likely not the first three propagating modes. 

The trouble with the results is mostly due to the following: a) there are errors 
in the discretization or b) the discretization is unable to properly model the boundary 
conditions. Indications of the first problem are evidenced in two ways. The first is the 
inability of the eigenvalues to converge to the correct analytical eigenvalues. Secondly, 
the plotting of the corresponding eigenvectors do not satisfy the boundary conditions 
as required. The solution to this problem is to find the correct way to difference the 
boundary conditions. Nothing can be done to fix the second problem, except to try 
a different discretization. The code as it stands does find some roots, but until the 
boundary conditions can be satisfied, there is no real confidence that the roots are 
comprehensive. 
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Several new problems arise from this research. The first is to determine a 
better way of differencing the boundary of the fluid-elastic layer. Once this problem 
is corrected, we can solve for the coefficients of transmitted and reflected waves given 
and incident wave. This is useful when extending the elastic case to the more complex 
porous case. If the porous case can be solved numerically, this will allow us to create 
ocean acoustic models where the transmitted and reflected waves can be determined 
from known coefficients. This has potential applications in sonar and mapping of 
the ocean floor. Finally, a similiar problem involves modeling a fluid layer only, with 
different sound velocity profiles and an acoustically hard boundary. The discretization 
would not involve the elastic boundary and may be a better way to approach the fluid- 
elastic case. 
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APPENDIX. MISCELLANEOUS DATA 



1. MATLAB CODE FOR GENERATING WAVE SPEEDS 

function [] =elastic(E, nu, rho) 

+ ************* + ******** ****** + * + * + *+ + 

'/, ELASTIC - This function computes the elastic constajits and 
% wavespeeds for various solids. Input values are the 

*/, Young's Modulus, Poisson's ratio, and the density. 

•/. 

y, Copyright 1997 by C. R. Myers, III 

lambda=0 ; 
mu=0; 
cl=0; 
ct=0 ; 

E=E+10“10; 

lambda = (E+nu)/((l+nu)*(l-2*nu)) 
mu = E/(2*(l+nu)) 

cl = round (sqrt( (lambda + 2*mu)/rho)) ; 
ct = round (sqrt (mu/rho)) ; 



2. DISCRETIZATION MATLAB CODE 

7o Fluid-Elasic Finite Difference Matlab Code 
7o 

®/o Copyright 1998 by Coley R. Myers, III 

% Initialize 
clear; 

format long e; 
w=0 ; 

E_val=[] ; 

'/, allows program to run in loop for obtaining 
'/, modes over a specified frequency range. 
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counter=l ; 



*/o loop is a counter and is used to obtain graph 
7 p graph of modes over a specified frequency range 

for loop=l : counter 

7p Vaxiables input for specific elastic solids 

I=sqrt(-1) ; 

Lf=10; 

Ls=200; 

w=2*pi+(w+. 1) ; rf=1026; rs=7700; 

ct=3145; cl=5690; cf=1500; 

modes=10; 

mu=rs*ct“2; 

lajn=rs*cl“2 - 2*mu; 

hf=l; 

hs=4; 

N=Lf /hf ; 

M=Ls/hs; 
d=N+M+M+2 ; 

A=zeros (d) ; 

kf=w/cf; kt=w/ct; kl=w/cl; 

y, Ensures matrices are clear when loop > 1 



Vl=[] ;Phi=[] ;Vphi=[] ;Psi=D ;Psixy=D ;Vphiy=[] ; 
V=[] ;D=[] ;Pres=[] ;U=D ; 
u=[] ;v=[] ;Styy=[] ;Stxx=D ;Stxy=D ; 
pr_styy=[] ;fl=D ;F1=D ;Y=D ;¥!=□ ; 

Vphia=[] ;Psia=D ;Phia=[] ; 

y, Fluid BC Phi(0)=0 and 1st equation 

p=i; 

A(p,l)=-(2/hf~2-kf"2) ; 

A(p,2)=l/hf'2; 

p=p+l; 

for i=2:N-l 
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A(p,i-l)=l/hf-2; 
A(p,i)=-(2/hf~2-kf~2) ; 
A(p, i+l)=l/hf “2; 
p=p+l ; 



end ; 

% 3rd BC for Phi 
A(p,N-l)=2/hf“2; 

A(p,N)=-2/hf“2+kf“2+(hs*w*2*rf)/(mu*hf) ; 
A(p,N+l)=(w*2/(hs*hf))*(l+hs*2*(kf'2-kl'2)) 
A (p ,N+2) =-hs*w“2*kl~2/hf ; 
A(p,N+3)=-w-2/(hs*hf ) ; 
A(p,N+M+2)=(kt“2*hs~2*w'2)/(2*hf) ; 
A(p,N+M+3)=-2*w~2/hf ; 

% 1st BC for Varphi 

p=p+l; 

A(p,N+M+2)=(2/hs)+kt~2*hs/2; 

A(p,N+M+3)=-2/hs; 

A(p,N+3)=-l/hs~2; 

A(p,N+2)=(2/hs'2)-kl~2; 

A(p,N+l)=kf'2-l/hs~2; 

A(p,N)=rf/mu; 

p=p+l; 

'/. 2nd BC for Varphi 

for i=(N+2):(N+M) 

A(p,i-l)=l/hs“2; 

A(p,i)=-(2/hs‘2-kl‘2) ; 

A(p,i+l)=l/hs'2; 

p=p+l; 

end; 

y, 3rd BC for Vcirphi 
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A(p,N+2*M+2)=(kf'2*hs/2)-2/hs; 

A(p,N+2*M+l)=2/hs; 

A(p,N+M+l)=kt‘2-l/hs‘2; 

A(p,N+M)=2/hs‘2-kl''2; 

A(p,N+M-l)=-l/hs“2; 

p=p+l; 



% 1st BC for Psi 

A(p ,N)=rf/ (mu*hs) ; 

A(p ,N+l)=kt~2/hs-2/hs''3; 
A(p ,N+2)=4/hs^3-2*kl^2/hs ; 
A(p,N+3)=-2/hs-'3; 
A(p,N+M+2)=2/hs''2; 
A(p,N+M+3)=-2/hs“2; 

p=p+l; 



y, 2nd BC for Psi 

for i=(N+M+3) : (M+2*M+1) 

A(p,i-l)=l/hs"2; 
A(p,i)=-(2/hs‘'2-kt“2) ; 
A(p,i+l)=l/hs"2; 
p=p+l; 

end; 

y. 3rd BC for Psi 

A(p,N+M-l)=2/hs“3; 
A(p,N+M)=(2+kl*2/hs)-4/hs^3; 
A(p,N+M+l) = (2/hs’'3)-kt"2/hs ; 
A (p ,N+2*M+1) =-2/hs“2 ; 
A(p,M+2*M+2)=2/hs~2; 

y. Final matrix is calculated 

[V, D]=eig(A); 
D=sqrt(diag(D)) ; 

[Y Il]=sort(D); 



% Lambda holds the number of modes specified for 
use with later calculations 

L 2 Lmbda=Y ( 1 : modes ) ; 

% when loop > 1, E_val holds values for graph 
% of modes over a specified frequency range 

E.val = [E.val Lajnbda] ; 

% Orders eigenvectors 

for i = Irmodes 
j=Il(i); 

V1=[V1 V(:J)] ; 

end 

y. Separates stress potentials for calculations 

for i = 1; modes 

Phi = [Phi Vl(l:N,i)] ; 

Vphi = [Vphi Vl(N+l:N+M+l,i)] ; 

Psi = [Psi Vl(N+M+2:d,i)] ; 
end; 

y. Solve for Pressure 

Pres = -rf.*Phi; 

U = -(l/(I*w)) .*Phi; 

y. Differencing for discretized potentials 

eta2 = Lambda. *2 - kl*2; 
xi2 = Lambda. *2 - kt*2; 
zeta2 = Lambda. *2 - kf"2; 

% Finding Psi_xy (same as Psi_yx) 

j=i; 

for i=l '.modes 
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Psixy(j ,i)=-((Psi(j+l,i)- (1 + (hs“2/2) ,*xi2(i)) .*Psi(j ,i)) ./ 
(hs + (hs“3/6) .*xi2(i))) ; 

end 

for j=2:M 

for i=l .-modes 

PsixyCj ,i)=(Psi(j-l,i)-Psi(j+l,i))/(2*hs) ; 

end 

end 



j=H+l; 

for i=l:modes 

Psixy(j,i) = (Psi(j-l,i)- (1 (hs~2/2).*xi2(i)).*Psi(j,i))./ 

(hs + (hs'‘3/6) .*xi2(i)) ; 

end 

"L Finding Vphi_y 

j=i; 

for i=l:modes 

VphiyCj ,i)=-((Vphi(j■^l,i)- (1 -t- (hs'‘2/2) . *eta2(i) ) . * 

Vphi(j ,i)) ./(hs + (hs“3/6) . *eta2(i) ) ) ; 

end 

for j=2:M 

for i=l :modes 

Vphiy ( j , i) = (Vphi ( j -1 , i) -Vphi ( j +1 , i) ) / (2*hs) ; 

end 

end 



for i=l:modes 

Vphiy (j , i)=(Vphi (j-1 , i)- (1 + (hs~2/2) .*eta2(i)) . ♦ 

Vphi (j , i) ) . /(hs (hs“3/6) .*eta2(i)) ; 

end 

*/. Solving for Tau_yy, Tau_xy, and Tau_yy 
for i=l:modes 

Styy = [Styy (-((lam + 2*mu) *kl“2) . *Vphi ( : , i) - (2*mu).* 
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Psixy(: ,i))] ; 

Stxx = [Stxx ((2*mu) . *Psixy( : , i) - ( (2*mu) *Lambda(i) *2 + 
lam*kl“2) .*Vphi(; ,i))] ; 

Stxy = [Stxy (((2*mu*I)*Lambda(i) .*Vphiy( : , i)) + 

(2*Lambda(i) "2 - kt“2)*(l/(I*Lambda(i))) .*Psi(: ,i))] ; 

end; 

*/, Solving for u and v 
for i=l:modes 

u = [u (I . *Lambda(i) . *Vphi ( : , i)+ (1 . / (I . *Lambda(i) ) ) . * 

PsixyC : ,i))] ; 

V = [v (Vphiy(:,i) - Psi(:,i))]; 

end 

% *********************** + ***:(:*****^**** + *****^:*** *♦*+:*** + + ***** 
y. Applying bi-orthogonality condition 

for i = Irmodes 

FI = [Pres(:,i); Styy(:,i)]; 
pr_styy = [pr.styy FI] ; 

end 

fl = (Stxx’*u - Stxy'*v) - Pres’*U; 

for j=l:modes 
for i=l:modes 

fl(i, j)=fl(i, j)*conj(fl(i,j)) ; 

end 

end 

y ♦*+***** + + + ** + + ** + ** + + + **+** + + ** + ****** ******** + =t:=M'*=t: + =l'*=t:**=(:** 
y Plot of norm of eigenvalues 

Yl=sqrt(Y. *conj (Y)) ; 

i=56:-l:l; 

figure 

plot(Yl(l:56),i, 'oO 
titleC'Norm of Eigenvalues 0 
xlabelC ’magnitude of lambda’) 
ylabel ( ’ modes ’ ) 
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^ 5^:5f:3jc3tc;t:5jc5jcst:5jcjf:3ic:jcsf:5^:;f:5jc5jc^3f:5jc3jc:4c3t::icsjc5ic:jc:4c*:ic:t:3tc5jc:jc:ic5tc:ic:jc3f:^:>tc^5)c;^:3jc3jcsf:st::f:sf:3tcjt:3f:3f:;t:5ic:ic3)c:*c;tc3tc3)c 

y. Determine coefficents for analytical solutions 



k=Lambda(l , 1) ; 

deta2 = k~2 - kl~2; 

dxi2 = k"2 - kt"2; 

dzeta2 = k"2 - kf"2; 

hs=200; 

hf=10; 

h=hf +hs ; 

y. Change value for C 
C=. 00004; 

Dl=(2*k''2-kt“2) ~2*tanh (sqrt (deta2) *(hs/2) ) ; 

D2=4*k“ 2*sqrt (deta2) *sqrt (dxi2) *t anh (sqrt (dxi2) * (hs/2) ) ; 

D3= (2=<=k~2-kt “2) "2*coth (sqrt (deta2) * (hs/2) ) =*=tajih (sqrt (dxi2) * 
(hs/2)); 

D= (D1-D2) / (D3-4*k~2*sqrt (deta2) *sqrt (dxi2) ) ; 
Bl=-(2+k"2-kt~2)/(2*I*k*sqrt (deta2)) ; 

B2= (cosh ( sqrt (dxi2)* (hs/2)) /cosh (sqrt(deta2)* (hs/2) ))*C; 
B=B1*B2; 

A2= (sinh (sqrt (dxi2) * (hs/2) ) / sinh (sqrt (deta2) * (hs/ 2) ) )*D ; 
A=B1*A2; 

E1=(D1-D2)*C; 

E2=ct/ (k*kt*sqrt (deta2) * (rf / rs) ) ; 

E3=co sh ( sqrt (dxi2)*hs/2) /sinh (sqrt (dzeta2)*(hf+hs/2)) ; 
E=E1*E2*E3; 

y=101:-l:l' ; 
yl=M+l:-l:l; 
y2=N:-l:l; 

Vphia = A . *cosh(sqrt (deta2. *yl) ) + B.*sinh(sqrt(deta2.*yl)) ; 
Psia = C.*cosh(sqrt(dxi2.*yl)) + D.*sinh(sqrt(dxi2.*yl)) ; 
Phia = E.*sinh(sqrt(dzeta2) .+(N-y2)) ; 

y, Plot analytical vs discretized solutions 
f igure 

Vphia=sqrt (Vphia. *conj (Vphia)) ; 
plot (Vphia, yl , 'o ^ ,Vphi( ; , 1) ,yl , 'xO 
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title( ’Analytical vs Discretized for Varphi (solid)’) 
xlabel(’x’) 

ylabel (’ analytical (o) - discretizes (x) ’ ) 

figure 

Psia = Psia.*I .*k; 

Psia = sqrt (Psia. *conj (Psia) ) ; 
plot (Psia, yl , ’o’ ,Psi(: ,1) ,yl, ’x’) 

title( ’Analytical vs Discretized for Psi (solid)’) 
xlabel ( ’ X ’ ) 

ylabel (’ ainalytical (o) - discretizes (x)’) 

figure 

Phia = I . *w. *Phia; 

Phia = sqrt (Phia. *conj (Phia) ) ; 
plot(Phia,y2,’o’ ,Phi ( : , 1) ,y2, ’x’ ) 

title( ’Analytical vs Discretized for Phi (fluid)’) 
xlabel(’x’) 

ylabel (’ analytical (o) - discretizes (x) ’ ) 

y, end loop 
end 
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